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ABSTRACT 

Aims. We aim to simulate radial profiles of molecular abundances and the gas temperature in cold and heavily shielded starless cores 
by combining chemical and radiative transfer models. Attention is also given to the time-evolution of both the molecular abundances 
and the gas temperature. 

Methods. A determination of the dust temperature in a modified Bonnor-Ebert sphere is used to calculate initial radial molecular 
abundance profiles. The abundances of selected cooling molecules corresponding to two different core ages are then extracted to 
determine the gas temperature at two time steps. The calculation is repeated in an iterative process yielding molecular abundances 
consistent with the gas temperature. Line emission profiles for selected substances are calculated using simulated abundance profiles. 
Results. The gas temperature is a function of time; the gas heats up as the core gets older because the cooling molecules are depleted 
onto grain surfaces. The change in gas temperature associated with depletion is of the order of 1 K. The contributions of the various 
cooling molecules to the total cooling power change with time, but the main cooling molecule at all times, in the range of environments 
studied here, is CO. Radial chemical abundance profiles are non-trivial: different species present varying degrees of depletion and in 
some cases inward-increasing abundances profiles, even at t > 10 5 years. Line emission simulations indicate that cores of different 
ages can present significantly different line emission profiles, depending on the tracer species considered. 

Conclusions. Chemical abundances and the associated line cooling power change as a function of time. Most chemical species are 
depleted onto grain surfaces at densities exceeding ~ 10 5 cm 4 . Notable exceptions are NH3 and N2IL; the latter is largely undepleted 
even at « H ~ 10 s cm" 3 . On the other hand, chemical abundances are not significantly developed in regions of low gas density even 
at t ~ 10 5 years, revealed by inward-increasing abundance gradients. Except in high-density regions where the gas-dust coupling is 
significant, the gas temperature can be significantly different from the dust temperature. This may have implications on core stability. 
Owing to the potentially large changes in line emission profiles induced by the evolving chemical abundance gradients, our models 
support the idea that observed line emission profiles can, to some extent, be used to constrain the ages of starless cores. 

Key words. ISM: abundances - ISM: clouds - ISM: molecules - Astrochemistry - Radiative transfer 



1. Introduction 

Low-mass stars are born out of gravitationally bound concen- 
trations of cold gas, the so-called starless cores. The very low 
temperature (~ 10K) of these objects allows molecular deple- 
tion onto grain surfaces to proceed efficiently, particularly in the 
dense central areas of the cores. Indeed, molecular depletion has 
been observed toward ma ny starless cores (e . g. | Willacv et alJ 
T998[ ICaselli et all fl999t iTafalla et all 1200a Uflreensen et alJ 
20041) . Because depletion is more efficient at higher densities, 



one expects the abundances of most molecules to fall toward the 
core center. Common molecules such as CO or NH3 may how- 
ever have (relatively) low abundances in young cores with low 
gas density, since in these conditions chemistry proceeds slowly. 
Therefore one can qualitatively expect the abundances of most 
species to peak at midrange (~ 10 3 - 10 4 cm" 3 ) densities, when 
studying regions where the gas density varies by multiple orders 
of magnitude. 

It has been observed and suggested by theoretical models 
that in starless cores the gas te mperature is generally different 
than t he dust temperature (e.g. Galli et al. 2002t lYoung et alJ 
12001 lKeto~& Field! 12005b iBergin et al] 120061) . This is because 
at low gas densities (corresponding to the outer parts of star- 
less cores) the gas-dust coupling is weak so that the gas cools 
efficiently through line radiation. Starless cores are supported 
(mainly) by thermal pressure; knowledge of the gas temperature 



is thus imperative to understanding the stability of the cores and 
hence the initial conditions of low-mass star formation. 

Because the line cooling power is determined by molecular 
abundances and the abundances themselves are time-dependent 
quantities, it is of interest to study how much the gas tempera- 
ture can change as a function of time and to what extent possi- 
ble temporal changes in gas temperature can feed back into the 
chemistry. Also, because observational data of starless cores is 
gathered through observations of line emission radiation, pre- 
dictions of radial molecular abundances, and the associated line 
emission, as a function of time are useful when interpreting ob- 
servations. In this paper, we combine chemical and radiative 
transfer calculations to study both the temporal and radial de- 
pendences of chemical abundances and the gas temperature in 
starless cores. 

The paper is structured as follows. In Sect.|2] we introduce 
the physical and chemical models and discuss the calculations 
and the various assumptions made. Section [3] presents the re- 
sults of our calculations. In Sect.|4]we discuss our results and in 
Sect.[5]we present our conclusions. 

2. Modeling 

In this section we give a detailed account of the physical and 
chemical models used. We also discuss the determination of the 
gas temperature and the spectral line simulations. 
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Table 1. Physical parameters of the MBESs considered in this 
paper. From left to right, the columns show the model denom- 
ination, core mass, core radius, total hydrogen density in the 
center shell of the model and the total hydrogen density in the 
outermost shell of the model, respectively. 



Model 


M[M G ] 


R oM [AU] 


« L H [cirr 3 ] 


7I° ut [CITT 3 ] 


A 


0.25 


3700 


1.85 x 10" 


2.19 x 10 5 


B 


1.0 


12600 


1.82 x 10 s 


1.95 x 10 4 


C 


5.0 


55000 


1.09 x 10 4 


1.02 x 10 3 



2. 1 . Physical model 

The physical core mo del is the modified Bonnor-Ebert 
sphere (hereafter M BES; iGalli et al.ll200l iKeto & Fieldll2005t 
ISipila et al.l |2011|) . The MBES is a gas sphere in hydro- 
static equilibrium bounde d by external p ressure. In contrast to 
the Bonnor-Ebert sphere (lBonnorlll956l) which is isothermal, 
the MBES presents a temperature struc ture that can be self- 
consistently calculated (ISipila et al.ll201 ll) . 

In this paper, we consider three different MBESs with masses 
0.25 Mq, I.OMq and 5.0 Mq. We chose to study several cores 
with different masses because of their varying physical proper- 
ties; low mass MBESs are dense and smal l, while high mas s 
MBESs are rather diffuse and more extended (Si pila et al.l 
From the chemical point of view this means that we expect to 
see very different radial profiles for species such as CO in the 
different cores, because molecular depletion is more efficient at 
higher densities. In all cases, we assume £1 = 6.0 for the non- 
dimensional r adius of the MBES , so that the model cores should 
be stable (see ISipila et alj|201 ll where the non-dimensional ra- 
dius £1 is defined). We assume the lOssenkopf & Hennindd 19941) 
thin ice mantle model for the optical properties of the dust and 
that the cores are embedded in a larger parent cloud with visual 
extinction corresponding to Ay = 10 at the core edge - with 
these choices, the model cores discus sed here correspond to the 
"type 1" cores of ISipila et al.l d201 ll) . The physical parameters 
of the model MBESs are presented in TableQ] the three different 
cores are denominated "A", "B" and "C". The (dust) temperature 
structures of the cores are shown in Fig.|2] 



2.2. Chemical model 

The chemical model program used in this paper utilizes the 
rate-equatio n method. The cod e is an updated version of the 
one used in ISipila et ail (1201 Ol) . expanded to include adsorp- 
tion, desorption and reactions on the surfaces of grains. The 
gas phase chemical reactions are adopted from the OSU re- 
action file osu_03_2QQ{fl. We assume spherical grains with 
fl g = 0.1 fim. The cosmic ray ionization rate is set to ( — 
1.3 x 10~ 17 s~'. The initial gas phase abundances and the sur- 
face reaction set (i.e. activation energ ies for select reactions) are 
adopted from ISemenov et al.l (120101) . The adopted initial abun- 
dances are atomic (with the exception of hydrogen which is prac- 
tically totally in H2) - the influence of this assumption on the 
results of this paper is discussed in Sect. 14.21 We assume that 
the total amount of matter on grain surfaces is initially zero, 
i.e. that the grains are bare. We note that this is a source of 
slig ht discrepancy with the dust t emperature calculation, since 
the lOssenkopf & Henningl (1 19941) dust model assumes grains 
coated in ice. 



The rate coefficient for adsorption of species i onto grains is 
&f ds = Vi crSi [cm 3 s ], where v,- = ^JW^T^Jmni is the thermal 
speed of species i, cr = na 2 g is the grain cross section and S , is 
a sticking coefficient, assumed here to be unity for all neutral 
species (we assume that ions do not stick). 

The ado pted desorption process is c osmic ray desorption as 
described in lHasegawa &Herbstl (1 19931) . i.e. the rate coefficient 
for cosmic ray desorption is kcs. = f ^td(70 K), where &td is 
the thermal desorption coefficient at 70 K and / is an efficiency 
factor (we assum e f _ 3.16 x 10~ 19 corresp onding to the value 
of £ adopted here; lHasegawa &~He rbst 1993). We do not include 
thermal desorption as a separate process because the core tem- 
peratures are low (< 10 K) - in these conditions, thermal desorp- 
tion is negligible compared to cosmic ray induced desorption. 
Also, the present model does not include photodesorption. 

The rate coefficients of reactions on grain surfaces are cal- 
culated inter nally in the chemical m odel program, following the 
formalism of lHasegawa et al.l (Il992l) . That is, the rate coefficient 
for a reaction between species i and j on the surface of a grain is 
calculated as 



kj = aK ij (R l m +R J m )/n i 



(1) 



where a is the branching ratio of the reaction and «d is the num- 
ber density of dust grains, k/j is an efficiency factor, assumed to 
be unity for exothermic reactions without activation energy. For 
endothermic reactions or exothermic reactions with activation 
energy, /c, ; - = exp(-£ , a /rd ust ) where £ a is the activation energy. 
The diffusion rate /?diff is calculated as 



Vo 

fidiff = — exp(-£ , diff /^ B rd) , 



(2) 



See http : //www . physics . ohio-state . edu/~eric/ 



where vo = ^2n s k^E< jTi 1 m is a characteristic vibration fre- 
quency, N s is the number of adsorption sites per grain and n s 
is the surface density of adsorption sites per grain. We assume 
n s = 1.5 X 10 15 cirr 2 , which leads to N s = 1.885 x 10 6 assuming 
spherical grains with a s =0.1 /mi. We assume E&g = 0.77 
for the diffus ion energies; mos t of th e binding energies E\, are 
adopted from lGarrod & Herbstl (120061) . but there are two excep- 
tions. ForH2 we have u sed a value E\, — 100 K; as pointed out by 
iGarrod & Paulvl (1201 ll) . the use of binding energies appropriate 
to a surface covered in water ice ca n lead to unphysical Hh ice 
abundances. IGarrod & Paulvl (1201 ll) corrected for this by allow- 
ing the binding energies of all species to change as the amount 
of H2 present on grain surfaces changes. Artificially lowering the 
H2 binding energy to Et, = 100K lowers the steady-state H2 ice 
abundance to about 10% of the H2O ice abundance, depending 
on the total hydrogen density. Our approach is rather arbitrary, 
but avoids the problem of increasing the number of equations to 
solve - decreasing the binding energy of H2 did not in our tests 
have a noticeable effect on the abundances of other species. For 
the bindi ng energy of atomic oxygen we have use d a value of 
1390 K dBergeron et alJl2008h ICazaux et alJl201 ll). Increasing 
the binding energy of O from 800 K (IGarrod & Herbstl 12006b 
to 1390 K changes greatly the behavior of O2 at high densi- 
ties; using the lower value leads to radial O2 abundances of 
~ 10~ 6 even in high density regions, inconsistent w ith observa- 
tions dGoldsmith et alJl2000t ILiseau. R. et al.ll2012l) . The higher 
binding energy value significantly decreases the radial O2 abun- 
dance at high densities (see Sect.[3] where the radial O2 abun- 
dances are further discussed). 

The effect of quantum tunneli ng on grai n surfa ces has 
been subject to some controversy. iKatz et all d 19991) argued 
that quantum tunneling of hydrogen on amorphous surfaces 
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would be inefficient. On the other hand, recent results of 
iGoumans & Anderssonl (1201 Ol) suggest that quantum tunneling 
dominates the CO + O — » CO2 reaction at low temperatures 
(this reaction is however very slow at typical starless core tem- 
peratures). In this work we do not include quantum tunnel- 
ing; this issue is further discussed in Sect. 14. II Also, we do not 
adopt the so called multilayer approach ( studied recently by e.g. 
iGarrod & Pauly|201 ltlTaquet et al. 2012) in which only the out- 
ermost surface layer is reactive. The exclusion of this process is 
also discussed in Sect.14.11 



2.3. Gas temperature calculation 

We have coupled radiative transfer calculations with a chemical 
model for a self-consistent determination of the gas tempera- 
ture ( our procedure is similar to the calculations of iBergin et al.l 
120061) . In practice, when the MBES has bee n constructed (using 
the method discussed in ISipila et ai]|2011l) . it is split into con- 
centric spherical shells by linearly dividing the total radius of 
the core to pi eces of equal leng th (the same approach was fol- 
lowed also in lSipila et al.ll201C4) . Consequently, each shell has a 
unique density and a unique temperature. Chemical evolution is 
then calculated separately in each shell, and the abundances of 
all substances as a function of time are extracted (initially, we 
assume r gas = r dust ). 

Using data from the chemistry model, molecular cooling cal- 
culations are performe d with a Monte Carlo radiative transfer 
program dJuveldfl997^ - the cooling calculations are similar 
to those presented in Juvela & Ysardl (1201 ll). i.e . we a lso em- 
ploy the gas-grain coupling scheme of lGoldsmithl (l200ll) . In this 
study, we consider the cooling rates of six substances: O, O2, C, 
12 CO, 13 CO and C ls O. CO is expected to be the main coolant 
- O and O2 are included to study their possible contribution to 
the total cooling power (see Sect. l3.2t . Isotopes are not explic- 
itly included in the chemical model - when inputting CO radial 
abundance profiles into the radiative transfer program, we have 
assumed 12 CO/ 13 CO and 12 CO/C 18 abundance ratios of 60 and 
500, respectiv ely (for typical values of the isotope ratios in the 
ISM, see e.g. IWilson & Roodlll994h . We also assume that the 
dust grains are not significantly cooled or heated in interactions 
with the gas, so that rd ust stays constant throughout the calcula- 
tion. 

Molecular abundances change with time - this means that 
the total cooling power may also change with time as molecules 
are depleted onto grain surfaces. To study how the gas temper- 
ature changes as a function of time owing to the changes in 
cooling power, we have run the cooling calculations using abun- 
dances corresponding to two different core ages (f = 10 5 and 
t = 10 years). The radial abundances as a function of time are 
further discussed in Sect.|3] 

After selecting the core age, the gas temperature as a func- 
tion of core radius is calculated. To assure consistency, the (first 
approximation) T g!is profile is then fed back into the chemi- 
cal model, which is now run assuming two separate temper- 
ature components. This yields new radial profiles for all sub- 
stances, which are input to another round of cooling calcula- 
tions. In practice only one or two iterations are needed - the 
cooling molecules are insensitive to small (~1-2K) changes in 
temperature and hence the gas temperature converges very fast 
toward a single solution. To illustrate this, we plot in the upper 
panel Fig. [T] the radial abundances of the cooling molecules at 
t = 10 6 years in model B as a function of core radius. Evidently 
the change in radial abundances is very small - after the ini- 
tial chemistry calculation (solid lines) the molecular abundances 




2000 4000 6000 8000 10000 12000 
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10 4 10 5 
t [years] 



Fig. 1. Upper panel: The abundances (with respect to «h, left- 
hand scale) of the cooling molecules (indicated in the figure) 
as a function of core radius in model B. Solid lines correspond 
to the initial chemistry calculation with r gas = 7d ust , other 
linestyles to subsequent iterations with r gas + T^x- The black 
curves represent r gas at different iterations (right-hand scale). 
Lower panel: The abundances of the cooling molecules and of 
OH as a function of time in a single-point chemical model with 



n H = 2.0 x 10 3 cm" 3 and T„ 



Tdu 



9K. 



change only at the edge of the core, and even there the change 
is negligible when further iterations are carried out (multiple 
curves with different linestyles practically superimposed). Also 
plotted in the figure are the gas temperatures at each iteration; 
the feedback from the minor changes in gas phase abundances 
to the gas temperature is negligible. 

As an example of the time-evolution of the cooling 
molecules, we plot in the lower panel of Fig. Q] the abundances 
of the cooling molecules in a single-point chemical model with 
total hydrogen density «h = 2 x 10 5 cm -3 and temperature 



?dust 



- gas 



9 K - these parameters represent roughly the cen- 



ter of model B. At this density, the CO abundance is decreasing 
at t — 10 s years due to depletion, but its abundance remains more 
or less constant from t — 10 s to t — 10 6 years. On the other hand, 
the abundances of C and O2 change by many orders of magni- 
tude due to depletion. The behavior of O2 around t — 10 5 years 
can be explained with the OH abundance. At earlier times, OH 
is primarily destroyed by N, C and O (in that order). After these 
atoms disappear from the gas phase, OH is destroyed mainly by 
H| and N2H + , but the destruction rate is an order of magnitude 
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lower. Both before and after t = 10 5 years, OH is mainly pro- 
duced in the dissociative electron recombination of H30 + ; the 
total production rate also decreases beyond 10 5 years, but the 
production rate decreases less than the destruction rate, leading 
to an OH abundance increasing with time. O2 is mainly produced 
in the neutral exchange reaction 

O + OH -> 2 + H (3) 
and mainly destroyed in 

C + 2 -> CO + O ; (4) 

because O disappears from the gas slightly faster than C, reaction 
© dominates over reaction © before 10 5 years and the abun- 
dance of O2 decreases. When the O abundance settles to a (more 
or less) constant value, the increasing OH abundance helps to 
increase also the O2 abundance. 

Figure Q] shows that the abundances of the cooling molecules 
can vary significantly with time. We shall see in Sect.[3]that this 
time evolution is reflected in the total cooling power of the core 
and hence on r gas . 

2.4. Simulated spectra 

Owing to the time evolution of the molecular abundances, one 
can expect line emission radiation observed toward starless cores 
to also change with time. That is, the line profiles of main tracer 
molecules such as N2H + and 13 CO are likely to be different de- 
pending on the core age. To quantify the effect of temporal vari- 
ation in molecular abundances on line emission radiation, we 
have calculated simulated line emission spectra using the ra- 
diative transfer program discussed above. We model the 1-0 
transitions of 13 CO, N 2 H + and HCO + as well as the (1,1) inver- 
sion transition of NH y Molecular collisio n data is taken from 
the LAMDA database dSchoier et alj|2005l) . 

In these calculations we study the model B core, assuming 
a distance of lOOpc. The beam FWHM is set equal to the core 
radius (126" at the assumed distance, cf. Table[TJ. We assume a 
turbulent linewidth of 0.22 km s~'. The results of these simula- 
tions are presented in Sect. 13.41 

3. Results 

3. 1 . Radial abundances of the cooling molecules 

Molecular abundances change with time, and this affects in par- 
ticular the total cooling power. We plot in the left-hand panels of 
Fig. [2] the radial abundances of the cooling molecules at differ- 
ent times in models A, B and C. It is evident that the abundances 
are strongly dependent on both time and density. For example, 
atomic carbon is more abundant than atomic oxygen in model A 
at t — 10 5 years, but the situation is reversed as the core grows 
older. At lower densities (models B and C), atomic carbon is gen- 
erally less abundant than atomic oxygen at both displayed times. 
The difference between the two species is particularly large in 
model C at t — 10 6 years. As illustrated by the lower panel of 
Fig.[T] the abundances of C and O can change drastically (ow- 
ing to depletion) in a relatively short time interval; this kind of 
behavior is accentuated in the low density of model C where 
chemical evolution is slow in general. Given enough time, the 
abundances of C and O would settle to similar levels even in 
model C. 

Molecular oxygen (green) evolves differently from atomic 
carbon and atomic oxygen; while its abundance decreases to- 
ward higher densities owing to depletion, the radial abundance 



of O2 is higher in an older core. At high densities, the growth 
of the radial abundance with time is about one order of mag- 
nitude, but up to two orders of magnitude at the low density of 
model C. It should be noted that at low density, the O2 abundance 
evolves differently than in the lower panel of Fig.[T| the local 
minimum present in Fig.[T]is, at low density, replaced by a local 
maximum. We note that our model C predicts sim ilar O2 abun- 
dances than the model of iHollenbach et all ( 120091) . even though 
they used the 800 K value for the binding energy of atomic oxy- 
gen. Furthermore, the O? abundances in model C agree with the 
abundances observed by iGoldsmith et al.l d201 lb , although they 
observed regions with much higher temperatures (~ 100 K) than 
in our models. 

CO (red) is fully depleted t — 10 s years in model A in the 
center of the core, but at the edge depletion is still underway. 
CO is fully depleted throughout the core by t — 10 6 years. In 
model B at t = 10 5 years, CO depletion has just begun in the 
center, but at the edge CO has not yet had time to reach its max- 
imum abundance (the evolution is similar to that shown in the 
lower panel of Fig.Q}. Again at t — 10 6 years CO is fully de- 
pleted throughout the core. The inward-increasing CO gradient 
seen in model C at t — 10 5 years is due to slower chemical evo- 
lution at the edge. In this model there is very little CO depletion 
at t — 10 6 years even in the core center. 

3.2. Gas temperatures and cooling rates 

The abundances of the cooling molecules vary considerably with 
time; the contributions of different molecules to the core cooling 
can change accordingly. Also, the total cooling power may de- 
crease as molecules are depleted onto grain surfaces. The right- 
hand panels of Fig.|2]plot the gas temperature in models A to C at 
two different times; evidently, the gas temperature changes with 
time. The changes in r gas as a function of time are the smallest 
in model A. The core is warmer at the edge at t — 10 6 years ow- 
ing to depletion; the main coolant CO is not yet fully depleted 
at the edge at t = 10 5 years. In models B and C the time vari- 
ation of r gas is more evident; however, in both cases the r gas 
time variation in a given point of the core is of the order of 1 K ; 
these results agree with the calculations of iBergin et al.l (120061) 
who also reported temporal r gas variations of ~1 K (albeit in a 
different scenario). 

To quantify the total cooling power and the contributions of 
the various molecules to it, we plot in Fig. [3] the cooling rates 
of the individual cooling molecules (indicated in the figure) and 
the total cooling rate as a function of gas density in model B at 
t = 10 5 and t = 10 6 years. There are clear differences between 
the two time steps. Atomic carbon is the second most impor- 
tant coolant at t — 10 5 years, but rather insignificant at t — 10 6 
years. The cooling rate of O2 is negligible and shows up only at 
the core edge at t = 10 6 years due to an increased abundance 
(cf. Fig.|2J. At this time, atomic carbon is significantly depleted 
which translates to a low cooling rate. Atomic oxygen is of little 
significance in terms of cooling power. 12 CO is the main coolant 
in both cases; its contribution totally dominates the total cooling 
power at later times. However, because of depletion, the abun- 
dance of CO drops by about an order of magnitude from t — 10 
to t — 10 6 years (cf. Fig.[2J. Accordingly, the total cooling rate 
is lower at the later time step and the core warms up as a result; 
a similar effect is seen in models A and C (Fig.|2j. 

lGoldsmifhl(l2001h argued that O and O2 would be insignif- 
icant coolants in dark clouds due to their low abundance. Our 
models confirm this statement. From our results it is also evi- 
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500 1000 1500 2000 2500 3000 3500 4000 ' 500 1000 1500 2000 2500 3000 3500 4000 



R [AU] R [AU] 




R [AU] R [AU] 

Fig. 2. Left-hand panels: Radial abundances (with respect to «h) of the cooling molecules (indicated in the figure) in models A 
(top), B (middle) and C (bottom). Solid lines correspond to t — 10 5 years, dashed lines to t — 10 6 years. Right-hand panels: The 
temperature profiles in models A to C. Solid lines correspond to t — 10 5 years, dashed lines to t = 10 6 years. Also plotted in each 
panel is the dust temperature (dotted lines). 

dent that the individual contributions of the various molecules to 3.3. Tracer molecules 
the core cooling evolve with time. 

Starless cores are often observed in various transitions of a mul- 
titude of molecules. This approach makes it possible to probe 
areas of different densities. However, as demonstrated by the 
cooling molecules above, the abundance of a given molecule 
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Fig. 3. The cooling rates of the cooling molecules (indicated in the figure) at t — 10 5 years (left panel) and at t — 10 6 years (right 
panel) in model B as a function of gas density. Also plotted in both panels is the total cooling power (black lines). 



can vary strongly with core radius, and most molecules can- 
not be used to probe the central parts of the cores because of 
depletion onto grain surfaces. Since the chemistry evolves rela- 
tively slowly in starless cores, we expect to see different abun- 
dance gradients, and particularly different degrees of depletion, 
in cores of different ages. 

To demonstrate this, we plot in Fig. [4] the radial abundance 
profiles of NH 3 , N 2 H + , CS, SO, HCO + and H 2 CO for two dif- 
ferent core ages in models A, B and C. It is evident that, in most 
cases, abundances change by about one order of magnitude when 
comparing the core center and edge, although the abundance gra- 
dients get somewhat less steep as the density increases. Notable 
exceptions are the sulfur-bearing molecules, the chemistry of 
which is however not very well understood due to a lack of ex- 
perimental data on rate coefficients. Core age is a very impor- 
tant factor. For example, in all models the abundances of NH3 
and N 2 H + are about an order of magnitude higher throughout 
the core at t — 10 6 years than at t — 10 5 years. HCO + displays 
similar behavior at low density, but at higher densities its radial 
abundance hardly changes with time. In contrast, the abundances 
of formaldehyde, SO and CS tend to drop as the cores get older. 
In models A and B, the sulfur-bearing molecules CS and SO are 
highly depleted at t = 10 6 years. 

Similarly to the cooling molecules, in model C the abun- 
dances are rather undeveloped even at t = 10 6 years due to 
slow chemical evolution. In this model the density is so low (cf. 
TableQ} that none of the species is significantly depleted even 
at t = 10 6 years. Indeed, model C is characterized by generally 
inward-increasing abundance gradients at t = 10 5 years. 

3.4. Line emission profiles 

The temporal variation of radial molecular abundances means 
from an observational point of view that observed line emission 
profiles depend on the core age, and hence line emission obser- 
vations can, at least in principle, be used in conjunction with the- 
oretical models to (roughly) estimate the age of a starless core. 

To demonstrate how the changes in radial abundances as a 
function of time translate to line emission radiation, we plot 
in Fig. [5] simulated line profiles of the 1-0 transitions of 



13 CO, N 2 H + and HCO + and the (1,1) inversion transition of 
NH3 in model B for two different core ages. The I3 CO pro- 
file is not significantly altered as the core grows older, although 
the peak intensity and optical depth decrease somewhat due to 
molecular depletion. The HCO + profiles present significant self- 
absorption; the optical depths at t — 10 5 and t = 10 6 years 
are 12.5 and 15.1, respectively. The increased optical depth at 
later times is due to a slightly larger abundance at the core edge 
(Fig.@). The N 2 H + and NH3 line profiles are significantly dif- 
ferent depending on assumed core age; this is again due to in- 
creased radial abundances (Fig.|4). 

Figure shows that there can be significant variation in line 
emission depending on the age of the studied core. For typical 
starless core densities ( 1 4 — 1 5 cm -3 ), NH3 and N 2 H + can serve 
as useful tracers for cores as old as t = 10 6 years. At even higher 
densities, represented here by model A, NH3 is depleted toward 
the core center but N 2 H + still retains a more or less constant 
radial abundance (Fig.|4}, which supports th e notion that N?H + 
is a viable tracer of high-density regions (e.g. lTafalla et alj2002t 
iPagani et al]l2007l) . 

Finally, we note that in the radiative transfer calculations, 
we have assumed zero infall velocity and that the core does not 
rotate. This means that the line profiles presented in Fig.|5]repre- 
sent an idealized case of a perfectly static core. As real cores are 
often either oscillat ing or contracting (e.g. lBroderick et al.l2008b 
lLee&Mversll2011h . care has to be taken when comparing simu- 
lated line profiles with observations. For example, the N 2 H + and 
NH3 profiles at t — 10 6 years shown in Fig.[5]are optically thick; 
the high optical thicknesses (8.12 and 7.63 for N 2 H + and NH3, 
respectively) are reduced to 6.06 and 6.61, respectively, if one 
assumes an infall velocity equal to the turbulent velocity used 
here (0.22 km s~'), producing sharper and slightly more intense 
line profiles. 

4. Discussion 

In this section, we discuss some aspects of chemical modeling 
that are not included in our calculations. We also discuss molec- 
ular abundances and the associated depletion predicted by our 
models. We end the section with a discussion on the issues in- 
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volved in determining gas temperatures and its possible impact 
on e.g. core stability. 

4.1. Quantum tunneling and multilayer grain chemistry 

We investigated the effect of quantum tunneling of atomic hydro- 
gen on gas phase abundances by running multiple single-point 
chemical models representing a range of densities, with H tun- 



neling either included o r excluded. In these test s, we set T^st = 
Tgas = 10 K. Following lHasegawa et alJ dl992h . whenever H is 
one of the reactants of a given surface reaction the k factor is 
calculated as 

k = exp[-2 (a /ti) (2 11 E c ,) 1/2 ] , (5) 

where a is an arbitrary barrier width and fi the reduced mass. 
lHasegawa et all (1 19921) assumed a = 1 A; we use this value 
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Fig. 6. The abundances (with respect to «h) of the cooling molecules (indicated in the figure) as a function of time in single-point 
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included; crosses: no tunneling. 
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as well. The value of a is a so urce of uncertain t y in c hemical 
models - recently, for example, iGarrod & Paulvl (1201 ll) used a 
value a — 2 A based on resul ts by lAndersson et al.l (120091) and 
iGoumans & Anderssonl(l2010t) . The diffusion rate of a tunneling 
reactant (in this case H) is calculated as 

*diff = W eX P I" 2 ™ C 2 m £ diff)' /2 ] " 

Figure [6] shows comparisons of the abundances of the cool- 
ing molecules as a function of time in a single-point chemical 
model with densities «h = 2 x 10 5 cm -3 (left-hand panel) and 
«h = 2 x 10 6 cm" 3 (right-hand panel). Evidently, the abundances 
of the cooling molecules are virtually unaffected by the inclusion 
of H tunneling. However, the inclusion of tunneling does affect, 
e.g., the abundances of hydrocarbons on grain surfaces - this 
is to be expected, since tunneling decreases the hydrogenation 
timescale on grain surfaces as hydrogen atoms scan the grain 
surfaces faster. A more extensive study of the effect of quantum 
tunneling is beyond the scope of this paper and it is unclear if in- 
cluding tunneling of, e.g., C and O would affect the abundances 
of the cooling molecules. We have performed the calculations 
presented in this paper without tunneling in order not to intro- 
duce many new unknown factors into the modeling. However, 
our tests indicate that at least in the case of H tunneling the re- 
sults of this paper should remain unaffected. Also, we do not 
expe ct including the so-called "modifi ed rate equation method" 
(e.g. lGarrodll2008t iDu & Parisell201 ]]) in the modeling of grain 
surface reactions to influence our results, as we consider rather 
large (0. 1 //m) grains in a low-temperature environment. 

As mentioned in Sect. 12.21 we do not consider in this pa- 
per the so-called multilayer approach, where only the outermost 
grain surface layer is reactive and the mantle under the su rface 
is chemically inert. As pointed out by iTaquet et al.l (1201 2|) , this 
means in particular that models such as the one considered here 
may underestimate the abundances of radicals on grain surfaces. 
It is unclear to what extent this affects gas phase abundances. 
There may be large differences between the multilayer and one- 
layer models at higher densities, where the depletion timescale is 
shorter. A quantitative study of this issue is certainly warranted, 
but beyond the scope of this paper. 

4.2. Molecular abundances and depletion 

We have constructed chemical models of starless cores assum- 
ing a modified Bonnor-Ebert sphere for the physical structure of 
the cores. All of the core models considered here present a den- 
sity gradient of about one order of magnitude; this property is 
imposed on the model cores due to the requirement of stability. 
To study larger density gradients, one could in principle consider 
supercritical core models with the non-dimensional critical ra- 
dius above £ ~ 6.5. However, in the context of supercritical and 
hence unstable cores it is probably not feasible to study model 
cores with ages much above ~ 10 5 years, since the core life- 
time is in this case expected to be very short compared to typical 
chemical timescales. 

In the core models considered here, we typically find molec- 
ular depletion of about one order of magnitude across the corefl 

2 This statement means that a given species is more depleted at the 
core center than at the edge by an order of magnitude. For example in 
model A at t = 10 6 years, CO is depleted by ~ 3 orders of magnitude 
at the core center and by ~ 2 orders of magnitude at the edge. Radial 
abundance profiles such as those presented here do not trace the "abso- 
lute" amount of depletion, but reflect the differential depletion arising 
from differences in density between different parts of the core. 



The depletion behavior of different species is however non- 
trivial: some species present flat radial abundance profiles while 
others are heavily depleted toward the core center. As demon- 
strated by Fig. [4] inward-increasing abundance gradients are pos- 
sible even in old (t ~ 10 6 years) cores, as long as the density is 
low. The density of the core plays an important role: abundance 
gradients are generally rather monotonous in high-density cores, 
while low-density cores can present more exotic abundance gra- 
dients because of the long timescales. 

In this study we have not attempted to identify viable tracer 
species; the abundance profiles presented in Sect. 13. 31 serve to 
demonstrate possible changes as a function of time that one 
could expect to observe in starless cores. In principle, models 
such as those presented here could be used for identifying tracer 
candidates, but this kind of analysis should be carried out using 
a wide parameter space, instead of fixing most physical param- 
eters such as grain size and cosmic ray ionization rate, as was 
done here. 

As noted in Sect. 12.21 in all models the gas phase abundances 
are initially atomic, with the exception of hydrogen which is al- 
most totally in H2. In this paper we have studied dense, deeply 
embedded cores - it is not clear to what extent the assumption of 
atomic initial abundances is applicable in this case. Starting from 
molecular instead of atomic abundances, i.e. assuming chemi- 
cal evolution previous to the starless core phase, will change the 
chemical evolution timescales. To study the possible effect of 
the initial abundances on the results of this paper, we have run 
multiple single-point chemical models assuming either atomic or 
molecular initial abundances. The molecular initial abundances 
were obtained by computing chemical evolution in a low density 
(~10 3 crrr 3 ) model; the so-obtained abundances were then used 
as input for higher density models. The atomic and molecular 
cases predict identical abundances after a time comparable to the 
age of the progenitor diffuse cloud - that is, if we let the diffuse 
model evolve for e.g. 10 6 years before extracting the abundances 
(i.e. initial abundances for the high density model), then the high 
density model predicts identical abundances at t > 10 6 years 
regardless of the assumed initial abundances. Thus, the results 
of this paper should hold to good accuracy whether the initial 
abundances are atomic or molecular, provided that the progen- 
itor diffuse core is not very old. We have refrained from more 
extensive testing of this issue because the present model does 
not allow for self-consistent merging of the diffuse and higher 
density models (by, e.g., allowing the diffuse cloud to gravita- 
tionally collapse into a denser configuration); certainly the state- 
ments made above should not be regarded as a general conclu- 
sion. For example in the context of collapse models, it has been 
show n that nitrogen che mistry is sensitive to the initial condi- 
tions dFlower et al.ll2006T). 

Recentlv. lHilv-Blant et all ( 1201 (j) and lPadovani et all (1201 ll) 
have studied observat ionally (and by using chemical models; 
iHilv-Blant et~ai1l2010h the abundances of CN, HCN and HNC 
in starless cores. They find that both HNC and HCN are present 
in appreciable abundances in the gas phase at densities where 
CO is depleted, and that the HNC/HCN ratio is close to unity. 
Our models reproduce both of these results. The HNC/HCN ra- 
tio depends on time; in young cores HNC is more abundant, but 
the situation is reversed as the core grows older. However, the 
ratio stays very close to unity at all times. For the abundances of 
CN, HCN and HNC we obtain values si milar (within an order o f 
magnitude) to the theoretical results of IHilv-Blant et all (1201 Oh . 
although they consider a collapse model whereas our model is 
static. I n conclusion, our model s als o fail to explain the ob serva- 
tions of lHilv-Blant et alJ (1201 Ol) and lPadovani et alJ d201ll) - the 
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possib le reasons for this failure are discussed in lHilv-Blant et all 
(l2010h . 



4.3. The NH 3 / N 2 H+ abundance ratio 

The NH3/N2H + abundance ratio has been studied in several 
cores and observed to increase towa rd the core centers (e.g. 
iTafalla et al.ll2002t iHotzel et al.ll2004l) . In these studies, the de- 
rived core densities are in the range 10 4 - 10 5 cirr 3 , corre- 
sponding roughly to our model B. Comparison between these 
observations and model B yields a rather good agreement at 
t = 10 5 years, our models reproducing rather well the radial 
abundances and consequently the NH3/N2H + abundance ratio, 
even though we do not attempt detailed modeling of any of 
the observed cores. Our models indicate that at high densities 
(> 10 5 cirT 3 ) the NH3/N2H + abundance ratio starts to fall as 
NH3 begins to deplete (while N2H + is still relatively unaffected). 
This is an interesting result and warrants observational verifica- 
tion. 

Similar modeling results of the NH3/N2H + ab undance ratio 
at inte rmediate densities have been published by lAikawa et al.l 
(120051) . In their models, however, N2H + depletes efficiently at 
higher densities, leading to a larger NH3/N2H + abundance ratio 
than in our models. The rea son for this discrepancy is proba- 
bly that lAikawa et aL I (120051) included in their reaction set the 
branching ratios of iGeppert et al.l (120041) for the N2H + + e~ re- 
combination reaction - these branching ratios have since been 
refuted and are not included in the osu_Q3_2008 reaction file. In 
the present model, the recombination yields N2 + H (90 %) and 
NH + N(10%). 



iFontani et al.l (1201 2l) have recently studied observationally 
the NH3/N2H + ratio in prestellar core candidates and found 
that the ratio increases strongly toward higher densities; in their 
Appendix A, they also show using a gas phase chemical model 
that the NH3/N2H + ratio should increase with increasing deple- 
tion. Our model seemingly predicts the opposite: the ratio de- 
creases at high density (corresponding to a higher degree of de- 
pletion). Howeve r, comparison of the results predicted by our 
model and that o f IFontani et al] d2012l) is difficult, because their 
model does not include a time-dependent treatment of freeze- 
out; instead, they set a value for the depletion factor fu describ- 
ing a uniform depletion of species heavier than helium. In our 
model, different species present highly variable degrees of de- 
pletion. This is demonstrated by e.g. the lower panel of Fig.Q] 
where CO is depleted by less than 2 orders of magnitude at 
t = 10 6 years, while e.g. C, O and N (not shown in the figure) are 
depleted by ~ 3 orders of magnitude. This differential depletion 
depends both on density and on the assumed core age - one can- 
not in our model assign a constant describing a general degree of 
depletion, and cons equently our r e sults c annot be directly com- 
pared with those of IFontani et al.l d2012l) . In our model, the ad- 
sorption of NH3 increases in importance as the density increases; 
this effect is primarily responsible for the decreasing NH3/N2H + 
ratio at high density. 

Finally, we no te that our chemica l model reproduces the 
modeling results of lFontani et al.l (1201 2l lower panel in their Fig. 
A-l) if we "turn off" the gas-grain interaction, i.e. if we consider 
gas phase chemistry and vary the metal abundances, although the 
dependence of the NH3/N2H + ratio on the depletion factor is in 
our model somewhat less steep than in IFontani et al.l d2012l) . 



4.4. Gas temperature 

The gas temperature is a function of time. As a core grows older, 
the total core cooling power decreases as cooling species (most 
importantly CO) are depleted onto grain surfaces and as a result, 
the gas temperature rises. However, in a typical core lifetime, 
this chan ge may be only ~ 1 - 2 K depending on the gas density 
(see also iBergin et alj|2006l) . Test calculations indicate that the 
feedback of small changes in gas temperature to the chemistry is 
small. It should be noted that this effect is not self-consistently 
modeled here as the gas temperature is not allowed to evolve in 
time; we have determined the gas temperature at different times 
using the gas phase abundances appropriate to those times. A 
more complete model taking into account dynamical changes in 
both the molecular abundances and the gas temperature could 
of course be built, but this approach is computationally much 
more demanding than the approach considered here and would 
probably not change our results significantly (as suggested by 
the rather small feedback of r gas on the chemistry). 

Gas temperature calculations are often carried out by ei- 
ther assuming constant radial abundances or assuming an ar- 
bitrary degree of depletion (using e.g. a step function) toward 
the core center that is the same for all cooling species (e.g. 
iKeto & Fieldl2005l:I.Tuvela & Ysardl201 lh . although calculations 
taking into account the chem istry have also been performed (e.g. 
IKeto & Casellill2008ll2010l) . Our results indicate that in reality 
the abundances are probably non-trivial; species deplete at dif- 
ferent rates and the degree of depletion varies from species to 
species. Furthermore, depletion tends to increase radially toward 
the core center. Overestimating the degree of depletion leads to 
a warmer core, while underestimating it leads to a colder core 
(see Sect. 13. 2\ . Figures|2]and|4]indicate that abundances in dense 
cores with advanced ages can be described by a (nearly) uniform 
degree of depletion, but this approximation breaks down for ei- 
ther young dense cores or cores with low gas densities. We have 
not made a direct comparison between our models and one as- 
suming e.g. constant radial profiles with arbitrary depletion near 
the center, but the differences in gas temperature predicted by the 
two approaches may not be larger than a few K - this of course 
also depends on the assumed cooling molecules. 

As a rule, the gas temperature in our models drops toward the 
core edge. This is because at the core edge the gas-dust coupling 
is weaker and t he photon escape prob ability higher than at the 
center (see also lJuvela & Ysardll201 ll) . We assume Ay = 10 at 
the edge, which means that heating of the cores by both the pho- 
toelectric effect and external UV radiation is negligible. In our 
models, there is no low-density gas outside the core that could 
provide heating by, e.g., diffusion. The presence of an external 
gas component would probably raise the outer core temperature 
somewhat; the chemistry would probably not be greatly affected 
(as discussed above), but a higher temperature could produce 
stronger line emission for the substances that are abundant in the 
outer, lower density regions (such as CO). 

In the calculations presented here, the physical core model 
is determined by the dust temperature and remains unchanged 
through subsequent chemistry and gas temperature calculations. 
This should be a good approximation in low mass (< 1 M Q ) 
cores, but at higher core masses the lower temperature of the 
gas in areas of moderate Ay (as opposed to that of the dust; 
see Fig.[2j may modify the density profile and affect the stabil- 
ity of the core (IBergin et al.ll2006t ISipila et al1l201ll) . although 
the presence of turbulent pressure (neglected here) could counter 
this. Furthermore, at low valu es of Ay, the gas heats up again due 
to photoelectric heating (e.g. IBergin et al.ll2006t IKeto & Casellil 
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120101) . providing additional support. Of course, changes in the 
density profile may also alter the chemistry, affecting the core 
cooling. Modeling these effects is beyond the scope of this pa- 
per, but a stud y of core stability including chemistry (expanding 
the analysis of Sipila et al] |201 ll) is planned. 

5. Conclusions 

We have calculated radial molecular abundances and gas temper- 
atures in core models representing starless cores in a semi- self- 
consistent way by inputting abundance profiles from a chemi- 
cal model into a radiative transfer program to determine the gas 
temperature. The gas temperatures were calculated at different 
timesteps corresponding to different stages of chemical evolu- 
tion. The gas temperature changes with time; the cores warm up 
slightly as cooling molecules are depleted onto grain surfaces. 
The interplay of evolving chemical abundances and the induced 
changes in gas temperature can translate into significant changes 
in line emission radiation as the core grows older. 

We examined the simulated radial abundances of some com- 
mon tracer molecules as a function of time in multiple core mod- 
els covering a range of densities. The models support earlier re- 
sults in the literature that NH3 and N2H + can resist depletion 
even at densities of 10 5 - 10 6 cnT 3 . On the other hand, many 
species present inward-increasing abundance gradients at low 
densities where the depletion timescale is long. In low density re- 
gions, chemical evolution is slow and significant depletion may 
not occur in a typical core lifetime. Molecular abundance gra- 
dients are in many cases non-trivial, especially for young cores 
with low densities. 

One interesting aspect of starless core chemistry, deutera- 
tion, is not covered in this paper. Deuteration is expected to 
be significant in the dense centers of the cores where heavy 
molecules are depleted onto grain surfaces and thus the deuter- 
ation process can proceed relatively unhindered. Because strong 
deuteration is effectively confined into a relatively small area in 
the core center, this process is not likely to have a notable in- 
fluence on the gas temperature and hence on the results of this 
paper. However, light deuterated molecules such as H2D + and 
D2H + are important tracers of high-density gas. Because of its 
importance in high-density regions, inclusion of deuteration to 
the chemistry model is underway. 
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